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Abstract 

Modeling stochasticity in gene regulatory networks is an important and complex problem in molecular systems 
biology. To elucidate intrinsic noise, several modeling strategies such as the Gillespie algorithm have been used 
successfully. This article contributes an approach as an alternative to these classical settings. Within the discrete 
paradigm, where genes, proteins, and other molecular components of gene regulatory networks are modeled as 
discrete variables and are assigned as logical rules describing their regulation through interactions with other 
components. Stochasticity is modeled at the biological function level under the assumption that even if the 
expression levels of the input nodes of an update rule guarantee activation or degradation there is a probability 
that the process will not occur due to stochastic effects. This approach allows a finer analysis of discrete models 
and provides a natural setup for cell population simulations to study cell-to-cell variability. We applied our methods 
to two of the most studied regulatory networks, the outcome of lambda phage infection of bacteria and the p53- 
mdm2 complex. 



1 Introduction 

Variability at the molecular level, defined as the phenoty- 
pic differences within a genetically identical population of 
cells exposed to the same environmental conditions, has 
been observed experimentally [1-4]. Understanding 
mechanisms that drive variability in molecular networks is 
an important goal of molecular systems biology, for which 
mathematical modeling can be very helpful. Different 
modeling strategies have been used for this purpose and, 
depending on the level of abstraction of the mathematical 
models, there are several ways to introduce stochasticity. 
Dynamic mathematical models can be broadly divided 
into two classes: continuous, such as systems of differential 
equations (and their stochastic variants) and discrete, such 
as Boolean networks and their generalizations (and their 
stochastic variants). This article will focus on stochasticity 
and discrete models. 

Discrete models do not require detailed information 
about kinetic rate constants and they tend to be more 
intuitive. In turn, they only provide qualitative informa- 
tion about the system. The most general setting is as fol- 
lows. Network nodes represent genes, proteins, and other 
molecular components of gene regulation, while network 
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edges describe biological interactions among network 
nodes that are given as logical rules representing their 
interactions. Time in this framework is implicit and pro- 
gresses in discrete steps. More formally, let X\, x n be 
variables, which can take values in finite sets X lt X n , 
respectively. Let X = Xi x ... x X n be the Cartesian pro- 
duct. A discrete dynamical system (DDS) in the variables 
a function 

f={fi fn):X^X 

where each coordinate function/: X — > X t is a function 
in a subset of {xi, x„). Dynamics is generated by itera- 
tion of /, and different update schemes can be used for 
this purpose. As an example, if X t = {0, 1} for all i, then 
each/ is a Boolean rule and/ is a Boolean network where 
all the variables are updated simultaneously. We will 
assume that each X t comes with a natural total ordering 
of its elements (corresponding to the concentration levels 
of the associated molecular species). Examples of this 
type of dynamical system representation are Boolean net- 
works, logical models and Petri nets [5-7]. 

To account for stochasticity in this setting several 
methods have been considered. Probabilistic Boolean 
networks (PBNs) [8,9] introduce stochasticity in the 
update functions, allowing a different update function to 
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be used at each iteration, chosen from a probability space 
of such functions for each network node. For other 
approaches, see [10-12]. These models will be discussed 
in more detail in the next section. In this article we pre- 
sent a model type related to PBNs, with additional fea- 
tures. We show that this model type is natural and a 
useful way to simulate gene regulation as a stochastic 
process, and is very useful to simulate experiments with 
cell populations. 

1.1 Modeling stochasticity in gene regulatory networks 

Gene regulation processes are inherently stochastic. 
Accurately modeling this stochasticity is a complex and 
important goal in molecular system biology. Depending 
on the level of knowledge of the biological system and 
the availability of data for it one could follow different 
approaches. For instance, viewing a gene regulatory net- 
work as a biochemical reaction network, the Gillespie 
algorithm can be applied to simulate each biochemical 
reaction separately generating a random walk corre- 
sponding to a solution of the chemical master equation 
of the system [13,14]. At an even more detailed level one 
could introduce time delays into the Gillespie simulations 
to account for realistic time delays in activation or degra- 
dation such as in circadian rhythms [15-17]. At a higher 
level of abstraction, stochastic differential equations [18] 
contain a deterministic approximation of the system and 
an additional random white noise term. However, all 
these schemes require that all the kinetic rate constants 
to be known which could represent a strong constraint 
due to the difficulty of measuring kinetic parameters, 
limiting these approaches to small systems. 

As mentioned in the introduction, discrete models are 
an alternative to continuous models, which do not 
depend on rate constants. In this setting, several 
approaches to introduce stochasticity have been pro- 
posed. Specially for Boolean networks, stochasticity has 
been introduced by flipping node states from 0 to 1 or 
vice versa with some flip probability [12,19-21]. However, 
it has been argued that this way of introducing stochasti- 
city into the system usually leads to over-representation 
of noise [11]. The main criticism of this approach is that 
it does not take into consideration the correlation 
between the expression values of input nodes and the 
probability of flipping the expression of a node due to 
noise. In fact, this approach models the stochasticity at a 
node regardless of the susceptibility to noise of the 
underlying biological function [11]. 

Probabilistic Boolean networks [8,9,22] is another sto- 
chastic method proposed within the discrete strategy. 
PBNs model the choice among alternate biological func- 
tions during the iteration process, rather than modeling 
the stochasticity of the function failure itself. We have 
adopted a special case of this setting, in which every 



node has associated to it two functions: the function 
that governs its evolution over time and the identity 
function. If the first is chosen, then the node is updated 
based on its logical rule. When the identity function is 
chosen, then the state of the node is not updated. The 
key difference to a PBN is the assignment of probabil- 
ities that govern which update is chosen. In our setting, 
each function gets assigned two probabilities. Precisely, 
let Xi be a variable. We assign to it a probability p: , 
which determines the likelihood that x t will be updated 
based on its logical rule, if this update leads to an 
increase/activation of the variable. Likewise, a probabil- 
ity pj determines this probability in case the variable is 

decreased/inhibited. The necessity for considering two 
different probabilities is that activation and degradation 
represent different biochemical processes and even if 
these two are encoded by the same function, their pro- 
pensities in general are different. This is very similar to 
what is considered in differential equations modeling, 
where, for instance, the kinetic rate parameters for acti- 
vation and for degradation/decay are, in principle, 
different. 

Note that all these approaches only take account of 
intrinsic noise which is generated from small fluctua- 
tions in concentration levels, small number of reactant 
molecules, and fast and slow reactions. Another source 
of stochasticity is related to extrinsic noise such as a 
noisy cellular environment and temperature. For more 
about intrinsic vs extrinsic noise see [3,23]. 

2 Method 

Our aim is to model stochasticity at the biological func- 
tion level under the main assumption that even if the 
expression levels of the input nodes of an update func- 
tion guarantee activation or degradation there is a prob- 
ability that the process will not occur due to 
stochasticity, for instance, if some of the chemical reac- 
tions encoded by the update function may fail to occur. 
This is similar to models based on the chemical master 
equation. This model type introduces activation and 
degradation propensities. More formally, let be 
variables which can take values in finite sets X lt X n , 
respectively. Let X = X\ x ... x X n be the Cartesian pro- 
duct. Thus, the formal definition of a stochastic discrete 
dynamical system (SDDS) in the variables 
collection of n triplets 

where 

♦ f t : X — > Xi is the update function for x b for all 
i = 1, n. 
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♦ p7 is the activation propensity. 

♦ pj is the degradation propensity. 

♦ p\> Pf 6 [0,1]- 

We now proceed to study the dynamics of such sys- 
tems and two specific models as illustration. 

2.1 Dynamics of SDDS 

Let F= {/i,/^,/^}" be a SDDS and consider x e X. 
For all i we define n it X (X[ — >fi{x)) and n ir x {x t — > xj) by 



2.1.1 Example 

Let n = 2, X = {0, 1} x {0, 1}, F = ifx,f 2 ): X -> X, where 
Table 4 represents the regulatory rules for xl and x2 
and Table 5 represents their propensity parameters 
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I- p], ifxi <fi{x), 

1 - pf, if Xj > 

1, if x, = /;(*)• 



That is, if the possible future value of the i-th coordi- 
nate is larger (smaller, resp.) than the current value, 
then the activation (degradation) propensity determines 
the probability that the i-th coordinate will increase 
(decrease) its current value. If the i-th coordinate and 
its possible future value are the same, then the i-th 
coordinate of the system will maintain its current value 
with probability 1. Notice that x (xi — > j,) = 0 for all 
Ji £ {xi,fi{x)}. 

The dynamics of F is given by the weighted graph X 
which has an edge from x e X to y e X if and only if 
Ji e {Xi, fi{x)} for all /. The weight of an edge x — > y is 
equal to the product 

n 

Wx-^y = Y\ 7T iA x i ->■ Yi) 

i=l 

By convention we omit edges with weight zero. See 
Additional file 1 for pseudocodes of algorithms to com- 
pute dynamics of SDDS. Software to test examples is 
available at http://dvd.vbi.vt.edu/adam.htmlt24] as a web 
tool (choose SDDS in the model type). 

Given F= j/i, p\,p\ J a SDDS, it is straightforward 

to verify that F has the same steady states (fixed points) 
as the deterministic system G = {/;}" =1 (see Additional 
file 1). It is also important to note that the dynamics of 
F includes the different trajectories that can be gener- 
ated from G using other common update mechanisms 
such as the synchronous and asynchronous schemes 
(see Additional file 1). 



Pr(01 -> 10) = (.1)(.9) = .09, Pr(01 -> 00) = (1 - 
.1)(.9) = .81 

Pr(01 -> 01) = (1 - .1)(1 - .9) = .09, Pr(01 -> 11) = 
(.1)(1 - .9) = .01 

Pr(10 -> 10) = (1 - .2)(1 - .5) = .4, Pr(10 -> 01) = 
(.2)(.5) = .1 

Pr(10 -> 00) = (.2)(1 - .5) = .1, Pr(10 -> 11) = (1 - .2) 
(.5) = .4 

Pr(ll -> 11) = (1)(1 - .9) = .1, Pr(ll 10) = (1) 
(.9) = .9 

Pr(00 H> 00) = (1)(1) = 1. 

Figure 1 shows that there is a 9% chance that the sys- 
tem will transition from 01 to 10. Similarly, there is an 
81% chance that the system will transition from 01 to 
00. The latter was expected because there is a high 
degradation propensity for f 2 . Note that 00 is a fixed 
point, i.e., there is 100% chance of staying at this state. 

3 Applications 

We illustrate the advantages of this model type by 
applying it to two widely studied biological systems, the 
regulation of the p53-mdm2 network and the control of 
the outcome of phage lambda infection of bacteria. 
These regulatory networks were selected because sto- 
chasticity plays a key role in their dynamics. 

3.1 Regulation in the p53-Mdm2 network 

The p53-Mdm2 network is one of the most widely studied 
gene regulatory networks. Abou-Jaude et al. [25] proposed 
a logical four-variable model to describe the dynamics of 
the tumor suppressor protein p53 and its negative regula- 
tor Mdm2 when DNA damage occurs. The wiring dia- 
gram of this model is represented in Figure 2, where P 
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Figure 1 State Space Diagram. This diagram depicts all trajectories to follow from any given initial state of the network. The numbers next to 
the edges specify the transition probabilities. Note that dynamics here is not deterministic, most of the states have different trajectories to 
follow. 



denotes cytoplasmic p53, nucleic p53, and the gene p53. 
Mc and Mn stand for cytoplasmic Mdm2 and nuclear 
Mdm2, respectively. DNA damage caused by ionic irradia- 
tion decreases the level of nucleic Mdm2 which enables 
p53 to accumulate and to remain active, playing a key role 
in reducing the effect of the damage. There is a negative 
feedback loop involving three components: p53 increases 
the level of cytoplasmic Mdm2 which, in turn, increases 
the level of nuclear Mdm2. Nucleic Mdm2 reduces p53 




Figure 2 Four-variable model for the p53-Mdm2 regulatory 
network. P, Mc, and Mn stand for protein p53, cytoplasmic Mdm2, 
and nuclear Mdm2, respectively. 



activity. This model also contains a positive feedback loop 
involving two components where p53 inhibits its negative 
regulator nucleic Mdm2. Note the dual role of P, as it 
positively regulates nucleic Mdm2 through cytoplasmic 
Mdm2. On the other hand, P negatively regulates nucleic 
Mdm2 by inhibiting Mdm2 nuclear translocation [25]. For 
more about the p53-Mdm2 system (see [4,25,26]). 

The dynamic behavior of the system is represented in a 
network of transitions called its state space (see Figure 3). 
This specifies the different paths to follow and the prob- 
abilities of following a specific trajectory from a given 
state. Dynamics here is not deterministic, i.e., most of the 
state vectors have different trajectories they can follow. 
The propensity parameters in Table 1 determine the likeli- 
hood of following certain paths. The state 0010 is a steady 
state, which is differentiated from the others by its oval 
shape. 

The state space for this model is specified by [0, 2] x 
[0, 1] x [0, 1] x [0, 1], that is, except for the first vari- 
able P which has three levels {0, 1, 2}, all other variables 
are Boolean. The update functions for this model are 
provided in Additional file 1 and also in the model 
repository of our web tool at http://dvd.vbi.vt.edu/adam. 
html. 

Individual cell simulations render plots similar to the 
ones shown in Figure 4. Each subfigure shows oscilla- 
tions as long as the damage is present with a variability 
in the timing of damage repair. On the other hand, cell 
population simulations, Figure 5, exhibit damped 
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Figure 3 State space diagram for parameters described in Table 1. The numbers next to the edges encode the transition probabilities. The 
order of the variables in each vector state is P, Mc, Mn, DNA damage. Self-loops are not depicted. States with darker background comprise the 
cycle with DNA damage. A second cycle with a lighter shaded background corresponds to the cycle with no DNA damage. The oval shaped 
state is a steady state. 



oscillations of the expression level of p53 as the degra- 
dation propensities of the damage increases. This is cor- 
related with the fact that, if the intensity of the damage 
is increased, more cells exhibit oscillations in the level 
of p53 which was experimentally observed in [4]. The 
initial state for all simulations was 0011 which repre- 
sents the state when DNA damage is introduced (0010 
is the steady state without perturbation). 

To highlight the features of our approach we compare 
our model with the one presented in [25] in which varia- 
bility has been analyzed. The main difference between 
these two models is in the way the simulations are 



Table 1 Propensity probabilities for the p53-Mdm2 
regulatory network 





P 


Mc 


Mn 


Dam 


Activation 


.9 


.9 


.9 


1 


Degradation 


.9 


.9 


.9 


.05 



Note that there is a low degradation propensity for DNA damage. 



performed. In [25], the transition from one state to the 
next is determined by parameters called "on" and "off" 
time delays. For instance, to transition from 2001 to 2101 
it is required that tjvic < tdam which means that the "on" 
delay for Mc (time for activating) is less than the "off 
delay (time for degrading) of the damage. Otherwise, if 
tMc > £ dam tri e system will transition from 2001 to 2000. 
In this article, transitions from one state to others are 
given as probabilities which are determined from the pro- 
pensity probabilities. Therefore, the complexity of the 
model presented here is at the level of the wiring diagram 
(i.e. the number of variables) while the complexity of the 
model in [25] is at the level of the state space (i.e. number 
of possible states) which is exponential in the number of 
variables. Another key difference is the way DNA damage 
repair is modeled. In [25], a delay parameter is asso- 
ciated with the disappearance of the damage, and this is 
decreased by a certain amount r at each iteration so that 

,(») _ M 
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Time steps Time steps 

Figure 4 Individual cell simulations for parameters described in Table 1 Each subfigure shows oscillations as long as the damage is 
present. This figure shows variability in the timing of damage repair and in the period of the oscilations. Each frame was generated from a 
single simulation with sixty time steps. The x-axis represents discrete time steps and the y-axis the expression level. The initial state for all 
simulations is 001 1 . 



iterations. In order to simulate DNA damage with this 

approach it is required to estimate r, n, and t^^. Within 

our model framework a single parameter, the degradation 
propensity, is used to model the damage repair which is a 
more natural setup. 

3.2 Phage lambda infection of bacteria 

Control of the outcome of phage lambda infection is 
one of the best understood regulatory systems [3,27,28]. 
Figure 6 depicts its core regulatory network that was 
first modeled by Thieffry and Thomas [28] using a logi- 
cal approach. This model encompasses the roles of the 
regulatory genes CI, CRO, CII, and N. From experimen- 
tal reports [3,28-30] it is known that, if the gene CI is 
fully expressed, all other genes are off. In the absence of 
CRO protein, CI is fully expressed (even in the absence 
of N and CII). CI is fully repressed provided that CRO 
is active and CII is absent. 

The dynamics of this network is a bistable switch 
between lysis and lysogeny, Figure 7. Lysis is the state 
where the phage will be replicated, killing the host. 
Otherwise, the network will transition to a state called 



lysogeny where the phage will incorporate its DNA into 
the bacterium and become dormant. It has been sug- 
gested [28,31] that these cell fate differences are due to 
spontaneous changes in the timing of individual bio- 
chemical reaction events. 

The state space for this model is specified by [0, 2] x 
[0, 3] x [0, 1] x [0, 1], that is, the first variable, CI, has 
three levels 0, 1, 2, the second variable, CRO, has four 
levels {0, 1, 2, 3}, and the third and fourth variables, CII 
and N, are Boolean. Update functions for this model are 
available in our supporting material, Additional file 1. 
This model has a steady state, 2000, and a 2-cycle invol- 
ving 0200 and 0300. The steady state 2000 represents 
lysogeny where CI is fully expressed while the other 
genes are off. The cycle between 0200 and 0300 repre- 
sents lysis where CRO is active and other genes are 
repressed. 

Cell population simulations were performed to mea- 
sure the cell-to-cell variability. Figure 8 was generated 
using the probabilities given in Tables 2 (top frame) and 
3 (bottom frame). The x-axis in both subfigures repre- 
sents discrete time steps while the y-axis captures the 
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Figure 5 Cell population simulations. Each subfigure was generated from 100 simulations, each representing a single cell with sixty time 
steps. Starting from the top left frame to the right bottom frame the degradation propensity for DNA damage was increased by 5%, i.e. 

Pdam = ' top ' ef=t '' Pdam = ' t0p rignt '' Pdam = (bottom left), and pj am = .2 (bottom right). The x-axis represents discrete 
time steps and the y-axis the average expression level. The initial state for all simulations was 001 1. This figure shows that, if the intensity of the 
damage is increased more cells exhibit oscillations in the level of p53, in agreement with experimental observations [4]. 



average expression level. The initial state for all simula- 
tions was 0000 which represents the state of the bacter- 
ium at the moment of phage infection. Figure 8 shows 
variability in developmental outcome, some of the net- 
works transition to lysis while others transition to lyso- 
geny. To measure how sensitive the dynamics of the 
network is to changes in the propensity probabilities, we 
have plotted the outcome of lysis-lysogeny percentages 
for different choices of these parameters. Figure 9 shows 
the variation in developmental outcome as a function of 
the propensity parameters of CI and CRO. Star points 
indicate the percentage of networks that transition to 
lysogeny and circle shaped points indicate the percentage 
of networks that end up in lysis. The bottom x-axis 



Table 2 Propensity parameters for Figure 7 (top frame) 





CI 


CRO 


CM 


W 


Activation 


.8 


.2 


.9 


.9 


Degradation 


.2 


.8 


.9 


.9 



There is a high activation propensity for CI while a low activation propensity 
for CRO. 



contains activation propensities for CI and degradation 
propensities for CRO while the top #-axis contains activa- 
tion propensities for CRO and degradation propensities 
for CI. The activation and degradation propensities for 
CU and N were all set equal to .9. Although the probabil- 
ity distributions for CI and CRO are very symmetric in 
Figure 9, it gives a good idea of how the variability in 
developmental outcome will change as the propensity 
parameters change. 

4 Conclusions 

Using a discrete modeling strategy, this article intro- 
duces a framework to simulate stochasticity in gene 



Table 3 Propensity parameters for Figure 7 (bottom 
frame) 





CI 


CRO 


CM 


W 


Activation 


.3 


.7 


.9 


.9 


Degradation 


.7 


.3 


.9 


.9 



There is a high activation propensity for CRO while a low activation 
propensity for CI. 
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regulatory networks at the function level, based on the 
general concept of PBNs. It accounts for intrinsic noise 
due to spontaneous differences in timing, small fluctua- 
tions in concentration levels, small numbers of reactant 



molecules, and fast and slow reactions. This framework 
was tested using two widely studied regulatory networks, 
the regulation of the p53-Mdrn2 network and the con- 
trol of phage lambda infection of bacteria. It is shown 
that in both of these examples the use of propensity 
probabilities for activation and degradation of network 
nodes provides a natural setup for cell population simu- 
lations to study cell-to-cell variability. The new features 
of this framework are the introduction of activation and 
degradation propensities that determine how fast or 
slow the discrete variables are being updated. This pro- 
vides the ability to generate more realistic simulations of 
both single cell and cell population dynamics. In the 
example of the p53-Mdm2 system, one can see that 
individual simulations show sustained oscillations when 
DNA damage is present, while at the cell population 
level these individual oscillations average to a damped 
oscillation. This agrees with experimental observations 
[4]. In the second example, A-phage infection of bac- 
teria, it is observed that differences in developmental 
outcome due to intrinsic noise can be captured with 
this framework. Due to the lack of experimental data we 
are unable to calibrate the model so that it reproduces 
the correct difference in percentages due to intrinsic 
noise. So instead we present a plot of the difference in 




Figure 7 State space for phage lambda model. The order of variables in each vector state is CI, CRO, Cll, N. The steady state 2000 represents 
lysogeny where CI is fully expressed while other genes are off. The cycle between 0200 and 0300 represents lysis where CRO is active and other 
genes are repressed. Self-loops are not depicted. 
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Figure 8 Cell population simulations. Both figures were generated from 100 simulations, each representing a single cell iteration often time 
steps. Top frame for parameters in Table 2 shows 93% lysis and 7% lysogeny while bottom frame for parameters in Table 3 shows 4% lysis and 
96% lysogeny. The x-axis represents discrete time steps while the y-axis shows the average expression level. The initial state for all the 
simulations is 0000. Solid (circle) points correspond to the average of CI (CftO), and dashed lines represent standard deviations. 



developmental outcome as a function of the propensity 
parameters. 

It is worth noting that this article addresses only 
intrinsic noise generated from small fluctuations in con- 
centration levels, small numbers of reactant molecules, 
and fast and slow reactions. Extrinsic noise is another 
source of stochasticity in gene regulation [3,23], and it 
would be interesting to see if this framework or a simi- 
lar setup can be adapted to account for extrinsic sto- 
chasticity under the discrete approach. This framework 
also lends itself to the study of intrinsic noise and it is 
useful for the study of developmental robustness. For 
instance, one could ask what the effect of this type of 
noise is on the dynamics of networks controlled by bio- 
logically inspired functions. 

Relating the propensity parameters to biologically mean- 
ingful information or having a systematic way for 



estimating them is very important. A preliminary analysis 
shows that it is possible to relate the propensity para- 
meters in this framework with the propensity functions in 
the Gillespie algorithm under some conditions (see Addi- 
tional file 1 where for a simple degradation model, the 
degradation propensity is correlated by a linear equation 
with the decay rate of the species being degraded). More 
precisely, in the Gillespie algorithm [13,14], if one discre- 
tizes the number of molecules of a chemical species into 
discrete expression levels such that within these levels the 
propensity functions for this species do not change signifi- 
cantly, then one obtains the setup of the framework pre- 
sented here as a discrete model. That is, simulation within 
the framework presented here can be viewed as a further 
discretization of the Gillespie algorithm, in a setting that 
does not require exact knowledge of model parameters. 
For a similar approach see [10]. 
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Figure 9 Variation in developmental outcome as a function of the propensity parameters Star points indicate the percentage of 
networks that transition to lysogeny and circle shaped points indicate the percentage of networks that end up in lysis. Bottom axis represents 
the activation (and degradation) propensities for CI (CRO) in increasing order. Likewise, the top axis represents the activation (and degradation) 
propensities for CRO (CI) in decreasing order. 



Acknowledgements 

DM and RL were partially supported by NSF grant 
CMMI-0908201. RL and DM thank Ilya Shmulevich for 
helpful suggestions. The authors thank the anonymous 
reviewers for many suggestions that improved the 
article. 

Additional material 



Additional file 1: Additional File 1 contains Supporting Material 



Author details 

'Department of Mathematics, Virginia Tech, Blacksburg, VA 24061-0123, USA 
2 Virginia Bioinformatics Institute, Virginia Tech, Blacksburg, VA 24061-0477, 
USA department of Mathematics, University of Nebraska-Lincoln, Lincoln, 
NE 68588, USA department of Computer Science, Virginia Tech, Blacksburg, 
VA 24061-0123, USA 



Competing interests 

The authors declare that they have no competing interests. 

Received: 14 November 2011 Accepted: 6 June 2012 
Published: 6 June 2012 

References 

1. E Avigdor, M Elowitz, Functional roles for noise in genetic circuits. Nature. 
467, 167-173 (2010). doi:l0.1038/nature09326 

2. M Acar, J Mettetal, A van Oudenaarden, Stochastic switching as a survival 
strategy in fluctuating environments. Nat Gen. 40, 471-475 (2008). 
doi:10.1038/ng.H0 

3. F St-Pierre, D Endy, Determination of cell fate selection during phage 
lambda infection. PNAS. 105, 20705-20710 (2008). doi:l0.1073/ 
pnas.0808831105 

4. N Geva-Zatorsky, N Rosenfeld, S Itzkovitz, R Milo, A Sigal, E Dekel, T 
Yarnitzky, Y Liron, P Polak, G Lahav, U Alon, Oscillations and variability in 
the p53 system. Mol Syst Biol. 2. 2006.0033, (2006) doi: 10.1038/msb4100068 

5. D Irons, Logical analysis of the budding yeast cell cycle. J Theor Biol. 257, 
543-559 (2009). doi:1 0.101 6/j.jtbi.2008.1 2.028 

6. R Thomas, R D'Ari, Biological Feedback (CRC Press, Boca Raton, 1990) 



Murrugarra ef al. EURASIP Journal on Bioinformatics and Systems Biology 2012, 2012:5 
http://bsb.eurasipjournals.eom/content/2012/1/5 



Page 11 of 1 1 



7. C Chaouiya, E Remy, B Mosse, D Thiery, Qualtative Aalysisof Regulatory 
Graphs: A Computational Tool Based on a Discrete Formal Framework, 
Lecture Notes in Control and Information Sciences, 294, 830-832 (2003) 

8. I Shmulevich, E Dougherty, S Kim, W Zhang, Probabilistic Boolean networks: 
a rule based uncertainty model for gene regulatory networks. 
Bioinformatics. 18(2), 261-274 (2002). doi:l 0.1 093/bioinformatics/l 8.2.261 

9. I Shmulevich, E Dougherty, Probabilistic Boolean Networks: The Modeling 
and Control of Gene Regulatory Networks, (SIAM, Philadelphia, 2010) 

10. S Teraguchi, Y Kumagai, A Vandenbon, S Akira, D Standley, Stochastic 
binary modeling of cells in continuous time as an alternative to 
biochemical reaction equations. Phys Rev E. 84(4) 062903 (2011) 

11. A Garg, K Mohanram, A Di Cara, G De Micheli, I Xenarios, Modeling 
stochasticity and robustness in gene regulatory networks. Bioinformatics. 
1525(12), il0l-i109 (2010) 

1 2. AS Ribeiro, SA Kauffman, Noisy attractors and ergodic sets in models of 
gene regulatory networks. J Theor Biol. 247, 743-755 (2007). doi:1 0.1 01 6/j. 
jtbi.2007.04.020 

13. D Gillespie, Exact stochastic simulation of coupled chemical reactions. J 
Phys Chem. 81(25), 2340-2361 (1977). doi:1 0.1 021 /j 1 00540a008 

14. D Gillespie, Stochastic simulation of chemical kinetics. Annu Rev Phys 
Chem. 58, 35-55 (2007). doi:l 0.1 146/annurev.physchem.58.032806.104637 

1 5. D Bratsun, D Volfson, LS Tsimring, J Hasty, delay-induced stochastic 
oscillations gene regulation. PNAS. 102(41), 14593-14598 (2005). 
doi:10.1073/pnas.0503858102 

16. AS Ribeiro, Stochastic and delayed stochastic models of gene expression 
and regulation. Math Biosci. 223(1), 1-11 (2010). doi: 1 0. 1 01 6/j. 
mbs.2009.1 0.007 

1 7. AS Ribeiro, R Zhu, SA Kauffman, A general modeling strategy for gene 
regulatory networks with stochastic dynamics. J Comput Biol. 13(9), 
1630-1639 (2006). doi:10.1089/cmb.2006.l3.1630 

18. T Toulouse, P Ao, I Shmulevich, S Kauffman, Noise in a small genetic circuit 
that undergoes bifurcation. Complexity. 11(1), 45-51 (2005). doi:10.1 002/ 
cplx.20099 

19. ER Alvarez-Buylla, A Chaos, M Aldana, M Benitez, Y Cortes-Poza, C Espinosa- 
Soto, DA Hartasanchez, RB Lotto, D Malkin, GJ Escalera Santos, P Padilla- 
Longoria, Floral morphogenesis: stochastic explorations of a gene network 
epigenetic landscape. PLoS ONF. 3(1 1), e3626 (2008). doi:l 0.1 371 /journal. 
pone.0003626 

20. Ml Davidich, S Bornholdt, Boolean network model predids cell cycle 
sequence of fission yeast. PLoS ONE. 3(2), e1672 (2008). doi:l 0.1 371 /journal. 
pone.0001672 

21. K Willadsen, J Wiles, Robustness and state-space structure of Boolean gene 
regulator models. J Theor Biol. 249(4), 749-765 (2008) 

22. R Layek, A Datta, R Pal, ER Dougherty, Adaptive intervention in probabilistic 
Boolean networks. Bioinformatics. 25(16), 2042-2048 (2009). doi:l0.1093/ 
bioinformatics/btp349 

23. SS Peter, BE Michael, DS Eric, Intrinsic and extrinsic contributions to 
stochasticity in gene expression. PNAS. 99(20), 12795-12800 (2002). 
doi:1 0.1 073/pnas.l 62041 399 

24. F Hinkelmann, M Brandon, B Guang, R McNeill, G Blekherman, A Veliz-Cuba, 
R Laubenbacher, ADAM: analysis of the dynamics of algebraic models of 
biological systems using computer algebra. BMC Bioinf 12, 295 (2011). 
doi:1 0.1 1 86/1471-21 05-1 2-295 

25. W Abou-Jaoude, D Ouattara, M Kaufman, From structure to dynamics: 
frequency tuning in the p53-mdm2 network: I. logical approach. J Theor 
Biol. 258(4), 561-577 (2009). doi:1 0.1 01 6/j.jtbi.2009.02.005 

26. E Batchelor, A Loewer, G Lahav, The ups and downs of p53: understanding 
protein dynamics in single cells. Nat Rev Cancer. 9, 371-377 (2009). 
doi:10.l038/nrc2604 

27. M Ptashne, A Genetic Switch: Phage X and Higher Organisms, (Cell Press 
and Blackwell Scientific Publications, Cambridge, 1992) 

28. D Thiery, R Thomas, Dynamical behaviour of biological regulatory networks- 
II. Immunity control in bacteriophage lambda. Bull Math Biol. 57, 277-295 
(1995) 

29. L Reichardt, D Kaiser, Control of A repressor synthesis. PNAS. 68, 2185-2189 
(1 971). doi:1 0.1 073/pnas.68.9.21 85 

30. P Kourilsky, Lysogenization by bacteriophage lambda I. Multiple infection 
and the lysogenic response. Mol Gen Gen. 122, 183-195 (1973). 

doi:1 0. 1 007/BF00435 1 90 



31. A Arkin, J Ross, H McAdams, Stochastic kinetic analysis of developmental 
pathway bifurcation in Phage X-infected Escherichia coll cells. Genetics. 149, 
1633-1648 (1998) 



doi:1 0.1 1 86/1 687-41 53-201 2-5 

Cite this article as: Murrugarra ef ah Modeling stochasticity and 
variability in gene regulatory networks. EURASIP Journal on Bioinformatics 
and Systems Biology 2012 2012:5. 



f \ 

Submit your manuscript to a SpringerOpen 0 
journal and benefit from: 

► Convenient online submission 

► Rigorous peer review 

► Immediate publication on acceptance 

► Open access: articles freely available online 

► High visibility within the field 

► Retaining the copyright to your article 



Submit your next manuscript at ► springeropen.com 

v ) 



